! This file calculates the tail area under the normal curve.
! Contributed by Alan Miller of CSIRO Division of Mathematics & Statistics, Clayton, Victoria.

module normcdf

implicit none

contains

	DOUBLE PRECISION FUNCTION NORMP(Z, P, Q, PDF)

!	Normal distribution probabilities accurate to 1.e-15.
!	Z = no. of standard deviations from the mean.
!	P, Q = probabilities to the left & right of Z.   P + Q = 1.
!       PDF = the probability density.
!       Based upon algorithm 5666 for the error function, from:
!       Hart, J.F. et al, 'Computer Approximations', Wiley 1968

!       Programmer: Alan Miller

!	Latest revision - 30 March 1986

	IMPLICIT DOUBLE PRECISION (A-H, O-Z)
	DATA P0, P1, P2, P3, P4, P5, P6 /220.2068679123761D0, &
      	  221.2135961699311D0, 112.0792914978709D0,       &
      	  33.91286607838300D0, 6.373962203531650D0,       &
      	  .7003830644436881D0, .3526249659989109D-01/
    DATA Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7 /440.4137358247522D0, &
      	  793.8265125199484D0, 637.3336333788311D0,           &
      	  296.5642487796737D0, 86.78073220294608D0,           &
      	  16.06417757920695D0, 1.755667163182642D0,           &
      	  .8838834764831844D-1/
    DATA CUTOFF /7.071D0/
	DATA ROOT2PI /2.506628274631001D0/

	ZABS = ABS(Z)

!	|Z| > 37.

	IF (ZABS .GT. 37.D0) THEN
	  PDF = 0.D0
	  IF (Z .GT. 0.D0) THEN
	    P = 1.D0
	    Q = 0.D0
		normp = P
	  ELSE
	    P = 0.D0
	    Q = 1.D0
		normp = P
	  END IF
	  RETURN
	END IF

!	|Z| <= 37.

	EXPNTL = EXP(-0.5D0*ZABS**2)
	PDF = EXPNTL/ROOT2PI

!	|Z| < CUTOFF = 10/sqrt(2).

	IF (ZABS .LT. CUTOFF) THEN
	  P = EXPNTL*((((((P6*ZABS + P5)*ZABS + P4)*ZABS + P3)*ZABS +   &
     		P2)*ZABS + P1)*ZABS + P0)/(((((((Q7*ZABS + Q6)*ZABS +   &
     		Q5)*ZABS + Q4)*ZABS + Q3)*ZABS + Q2)*ZABS + Q1)*ZABS +  &
     		Q0)
	  normp = P

!	|Z| >= CUTOFF.

	ELSE
	  P = PDF/(ZABS + 1.D0/(ZABS + 2.D0/(ZABS + 3.D0/(ZABS + 4.D0/  &
     		(ZABS + 0.65D0)))))
	  normp = P
	END IF

	IF (Z .LT. 0.D0) THEN
	  Q = 1.D0 - P
	ELSE
	  Q = P
	  P = 1.D0 - Q
	  normp = P
	END IF
	RETURN
	END FUNCTION NORMP
	
	! ----------------------------------------------------------------------
	! The routines for random number generation were obtained from
	! the book NUMERICAL RECIPES


	SUBROUTINE ran1sub(idum,x) 
	IMPLICIT NONE
	INTEGER, PARAMETER :: K4B=selected_int_kind(9)
	INTEGER(K4B), INTENT(INOUT) :: idum
	REAL(8) :: x 
	INTEGER(K4B), PARAMETER :: IA=16807, IM=2147483647,IQ=127773,IR=2836
	REAL(8), SAVE :: am
	INTEGER(K4B), SAVE :: ix=-1,iy=-1,k
	if (idum <= 0 .or. iy < 0) then
		 am=nearest(1.0,-1.0)/IM
		 iy=ior(ieor(888889999,abs(idum)),1)
		 ix=ieor(777755555,abs(idum))
		 idum=abs(idum)+1
	end if
	ix=ieor(ix,ishft(ix,13))
	ix=ieor(ix,ishft(ix,-17))
	ix=ieor(ix,ishft(ix,5))
	k=iy/IQ
	iy=IA*(iy-k*IQ)-IR*k
	if (iy<0) iy=iy+IM
	x=am*ior(iand(IM,ieor(ix,iy)),1)
	END SUBROUTINE ran1sub


	SUBROUTINE ran1(idumtemp,vv,n)

	implicit none
	!external ran1sub

	REAL(8) vv(n),x
	INTEGER idumtemp,kk,n

	vvloop: do kk = 1,n
		 call ran1sub(idumtemp,x)
		 vv(kk) = x
	end do vvloop

	END SUBROUTINE ran1


	!modified from http://jean-pierre.moreau.pagesperso-orange.fr/Fortran/tmdian_f90.txt May 3, 2015
	!***************************************************************
	!*  Calculate the median value of an array with the Heapsort   *
	!*  method                                                     *
	!* ----------------------------------------------------------- *
	!* REFERENCE:                                                  *
	!*      "NUMERICAL RECIPES By W.H. Press, B.P. Flannery,       *
	!*       S.A. Teukolsky and W.T. Vetterling, Cambridge         *
	!*       University Press, 1986" [BIBLI 08].                   *
	!* ----------------------------------------------------------- *
	!*******************************************************
	!* Given an array X of N numbers, returns their median *
	!* value XMED. The array X is modified and returned    *
	!* sorted in ascending order.                          *
	!*******************************************************
	SUBROUTINE PERCENTILES(X,N,YP10,YP25,YMED,YP75,YP90)
	  implicit none
	  integer N, N10, N2, N4, N34, N90
	  double precision X(N), YMED, YP10, YP25, YP75, YP90
	  call hpsort(N,X)
	  N10 = N/10
	  N2=N/2
	  N4=N/4
	  N34=3*N/4
	  N90 = 9*N/10
	  !median
	  if (2*N2.eq.N) then
	    YMED = 0.5*(X(N2)+X(N2+1))
	  else
	    YMED = X(N2+1)
	  endif
	  !p10 - rough
	  YP10=X(N10+1)
	  !p25 - rough
	  YP25=X(N4+1)
	  !p75 - rough
	  YP75=X(N34+1)
	  !p90 - rough
	  YP90=X(N90+1)
	  return
	END SUBROUTINE PERCENTILES


	!*****************************************************
	!*  Sorts an array RA of length N in ascending order *
	!*                by the Heapsort method             *
	!* ------------------------------------------------- *
	!* INPUTS:                                           *
	!*	    N	  size of table RA                       *
	!*      RA	  table to be sorted                     *
	!* OUTPUT:                                           *
	!*	    RA    table sorted in ascending order        *
	!*                                                   *
	!* NOTE: The Heapsort method is a N Log2 N routine,  *
	!*       and can be used for very large arrays.      *
	!*****************************************************         
	SUBROUTINE HPSORT(N,RA)
	  implicit none
	  integer N, L, IR, J, I
	  double precision RA(N), RRA
	  L=N/2+1
	  IR=N
	  !The index L will be decremented from its initial value during the
	  !"hiring" (heap creation) phase. Once it reaches 1, the index IR 
	  !will be decremented from its initial value down to 1 during the
	  !"retirement-and-promotion" (heap selection) phase.
	10 continue
	  if(L > 1)then
	    L=L-1
	    RRA=RA(L)
	  else
	    RRA=RA(IR)
	    RA(IR)=RA(1)
	    IR=IR-1
	    if(IR.eq.1)then
	      RA(1)=RRA
	      return
	    end if
	  end if
	  I=L
	  J=L+L
	20 if(J.le.IR)then
	  if(J < IR)then
	    if(RA(J) < RA(J+1))  J=J+1
	  end if
	  if(RRA < RA(J))then
	    RA(I)=RA(J)
	    I=J; J=J+J
	  else
	    J=IR+1
	  end if
	  goto 20
	  end if
	  RA(I)=RRA
	  goto 10
	END SUBROUTINE HPSORT
	


!Subroutine to find the inverse of a square matrix
!Author : Louisda16th a.k.a Ashwith J. Rego
!Reference : Algorithm has been well explained in:
!http://math.uww.edu/~mcfarlat/inverse.htm           
!http://www.tutor.ms.unimelb.edu.au/matrix/matrix_inverse.html
SUBROUTINE FINDInv(matrix, inverse, n, errorflag)
	IMPLICIT NONE
	!Declarations
	INTEGER, INTENT(IN) :: n
	INTEGER, INTENT(OUT) :: errorflag  !Return error status. -1 for error, 0 for normal
	DOUBLE PRECISION, INTENT(IN), DIMENSION(n,n) :: matrix  !Input matrix
	DOUBLE PRECISION, INTENT(OUT), DIMENSION(n,n) :: inverse !Inverted matrix
	
	LOGICAL :: FLAG = .TRUE.
	INTEGER :: i, j, k, l
	REAL :: m
	REAL, DIMENSION(n,2*n) :: augmatrix !augmented matrix
	
	!Augment input matrix with an identity matrix
	DO i = 1, n
		DO j = 1, 2*n
			IF (j <= n ) THEN
				augmatrix(i,j) = matrix(i,j)
			ELSE IF ((i+n) == j) THEN
				augmatrix(i,j) = 1
			Else
				augmatrix(i,j) = 0
			ENDIF
		END DO
	END DO
	
	!Reduce augmented matrix to upper traingular form
	DO k =1, n-1
		IF (augmatrix(k,k) == 0) THEN
			FLAG = .FALSE.
			DO i = k+1, n
				IF (augmatrix(i,k) /= 0) THEN
					DO j = 1,2*n
						augmatrix(k,j) = augmatrix(k,j)+augmatrix(i,j)
					END DO
					FLAG = .TRUE.
					EXIT
				ENDIF
				IF (FLAG .EQV. .FALSE.) THEN
					PRINT*, "Matrix is non - invertible"
					inverse = 0
					errorflag = -1
					return
				ENDIF
			END DO
		ENDIF
		DO j = k+1, n			
			m = augmatrix(j,k)/augmatrix(k,k)
			DO i = k, 2*n
				augmatrix(j,i) = augmatrix(j,i) - m*augmatrix(k,i)
			END DO
		END DO
	END DO
	
	!Test for invertibility
	DO i = 1, n
		IF (augmatrix(i,i) == 0) THEN
			PRINT*, "Matrix is non - invertible"
			inverse = 0
			errorflag = -1
			return
		ENDIF
	END DO
	
	!Make diagonal elements as 1
	DO i = 1 , n
		m = augmatrix(i,i)
		DO j = i , (2 * n)				
			   augmatrix(i,j) = (augmatrix(i,j) / m)
		END DO
	END DO
	
	!Reduced right side half of augmented matrix to identity matrix
	DO k = n-1, 1, -1
		DO i =1, k
		m = augmatrix(i,k+1)
			DO j = k, (2*n)
				augmatrix(i,j) = augmatrix(i,j) -augmatrix(k+1,j) * m
			END DO
		END DO
	END DO				
	
	!store answer
	DO i =1, n
		DO j = 1, n
			inverse(i,j) = augmatrix(i,j+n)
		END DO
	END DO
	errorflag = 0
END SUBROUTINE FINDinv


	
	
end module normcdf